-------------------------------------------------------------------------------
---
-- Example of use of the Ada Binding to MPFR, the Multiple Precision Floating
-- Point Reliable Library. See <http://www.mpfr.org/>.
--
-- Copyright (C) 2009 Vincent DIEMUNSCH.
-- GNU General Public License.
--
-- This program is free software: you can redistribute it and/or modify it under
-- the terms of the GNU General Public License as published by the Free Software
-- Foundation, either version 3 of the License, or (at your option) any later
-- version.
--
-- This file is distributed in the hope that it will be useful, but WITHOUT ANY
-- WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-- A PARTICULAR PURPOSE. See the GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License along with
-- this program.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------------------
---
--
-- This little program is an example of computing an approximation of
-- 1+1/1!+1/2!+...+1/100! using the "thick" Ada binging to MPFR: MPFR.Floats
-- To make a comparison with the "thin" Ada binding, MPFR, we give also
-- the "translation" in Ada of the sample C program given by the MPFR site,
-- see : < http://www.mpfr.org/sample.html >
--
-- The result of this program is (for 61 digits, and 204 bits) :
--    Sum is 2.718281828459045235360287471352662497757247093699959574966968
--

with Ada.Text_IO;      use Ada.Text_IO;
with Interfaces.C;     use Interfaces.C;
with GMP.Integers;
with MPFR.Floats;
with Rational_Numbers; use Rational_Numbers;

procedure MPFR_Example is

   Precision : constant := 61; -- Digits.

begin
   New_Line;
   Put_Line ("Example of use of the Ada Binding to MPFR, the Multiple Precision "
            &"Floating Point Reliable Library.");
   Put_Line ("Copyright (C) 2009 Vincent DIEMUNSCH.");
   Put_Line ("GNU General Public License.");
   New_Line;

   declare
      R : Rational;
   begin
      R := 1 + 5 / (-3);
      Put_Line (" R := " & Image (R));
      Put_Line (" trunc := " & Integer'Image (Rational_Numbers.int (1 + 5 / (-3))));
      Put_Line (" floor := " & Integer'Image (Rational_Numbers.Floor (1 + 5 / (-3))));
   end;


   Put_Line ("Computing sum of 1/n! for n in 0 .. 100");
   New_Line;

   -- 1) Thick binding version :
   -- Note that we use Integers because they are EXACT values, to initialize
   -- Multiple precision floats or in computations with the latters
   -- and not float values, that are approximations of real values with a too
   -- small precision.
   declare
      use MPFR.Floats;
      Fact, Sum : MPFR_Float;
   begin
      Put_Line ("1. Use of MPFR_Floats with a precision of : " &
               Integer'Image (Precision) & " digits.");
      Set_Default_Precision (Precision);
      Put_Line ("   Default precision is set to            : " &
                Integer'Image (Get_Default_Precision) & " digits.");
      Fact := To_MPFR_Float (1);
      Sum  := To_MPFR_Float (1);
      for i in 1 .. 100 loop
         Fact := i * Fact;
         Sum  := Sum + 1 / Fact;
      end loop;
      Put_Line ("Sum is" & Image (Sum));
      New_Line;
   end;

   -- 2) Thin binding version : a simple "translation" of the C program,
   --    except the fact that we use here Round to Nearest, instead of
   --    round to minus infinity. We take the default precision set by
   --    the Thick Binding version.
   declare
      use MPFR;
      s, t, u : mpfr_t;
   begin
      Put_Line ("2. Direct call of MPFR with a precision of : " &
                 mpfr_prec_t'Image (MPFR.mpfr_get_default_prec) & " bits.");
      mpfr_init (t);
      mpfr_set_d (t, 1.0, GMP_RNDN);
      mpfr_init (s);
      mpfr_set_d (s, 1.0, GMP_RNDN);
      mpfr_init (u);
      for i in 1 .. 100 loop
         mpfr_mul_ui (t, t, unsigned_long (i), GMP_RNDN);
         mpfr_set_d (u, 1.0, GMP_RNDN);
         mpfr_div (u, u, t, GMP_RNDN);
         mpfr_add (s, s, u, GMP_RNDN);
      end loop;
      Put ("Sum is ");
      mpfr_out_str (null, 10, Precision, s, GMP_RNDN);
      New_Line;
      mpfr_clear (s);
      mpfr_clear (t);
      mpfr_clear (u);
      New_Line;
   end;

   -- 3) Check with floats, Long_Floats and Long_Long_Floats
   declare
      Fact, Sum                     : Float;
      Long_Fact, Long_Sum           : Long_Float;
      Long_Long_Fact, Long_Long_Sum : Long_Long_Float;
      use MPFR.Floats;
      MPFRSum                       : MPFR_Float;
   begin
      Put ("3. Check result with floats :               ");
      Fact := 1.0;
      Sum := 1.0;
      for i in 1 .. 100 loop
         Fact := Float (i) * Fact;
         Sum  := Sum + 1.0 / Fact;
      end loop;
      Put_Line (Float'Image (Sum));

      Put ("   Then check with Long_Floats :            ");
      Long_Fact := 1.0;
      Long_Sum := 1.0;
      for i in 1 .. 100 loop
         Long_Fact := Long_Float (i) * Long_Fact;
         Long_Sum  := Long_Sum + 1.0 / Long_Fact;
      end loop;
      Put_Line (Long_Float'Image (Long_Sum));

      Put ("   And Finaly check with Long_Long_Floats : ");
      Long_Long_Fact := 1.0;
      Long_Long_Sum := 1.0;
      for i in 1 .. 100 loop
         Long_Long_Fact := Long_Long_Float (i) * Long_Long_Fact;
         Long_Long_Sum  := Long_Long_Sum + 1.0 / Long_Long_Fact;
      end loop;
      Put_Line (Long_Long_Float'Image (Long_Long_Sum));

      MPFRSum        := To_MPFR_Float (Sum);
      Put_Line
        ("   Float value converted to a MPFR :        " & Image (MPFRSum));
   end;

   New_Line;

   -- 4) compare precision on the Golden Number :
   declare
      use MPFR.Floats;
      Golden_Long_Float : Long_Float := 1.0;
      Golden_MPFR_Float : MPFR_Float := To_MPFR_Float (1);
   begin
      for i in 1 .. 100 loop
         Golden_Long_Float := 1.0 + 1.0 / Golden_Long_Float;
         Golden_MPFR_Float := To_MPFR_Float (1) +
                              To_MPFR_Float (1) / Golden_MPFR_Float;
      end loop;
      Put_Line
        ("Golden Number computed  in Long Float: " &
         Long_Float'Image (Golden_Long_Float));
      Put_Line
        ("Golden Number converted to Long Float: " &
         Long_Float'Image (Round(Golden_MPFR_Float)));
      Put_Line
        ("Golden Number in Multiprecision Float: " & Image (Golden_MPFR_Float));
      Put_Line
        ("Error                                : " &
         Long_Float'Image (Compare (Golden_Long_Float, Golden_MPFR_Float)));
   end;

   -- 5) With Unbounded Integers
   declare
      use GMP.Integers;
      N : Unbounded_Integer := To_Unbounded_Integer(
      "12345678901234567890123456789012345678901234567890123456789012345678901234567890123456");
      use MPFR.Floats;
      X : MPFR_Float := To_MPFR_Float (N, Best_Precision);

   begin
      New_Line;
      Put_Line ("N =  "& Image (N));
      Put_Line ("X ="  & Image (X));
   end;

   -- 6) Test reported by stef.van-vlierberghe@eurocontrol.int, Oct 5, 2011

--     declare
--        use MPFR.Floats;
--     begin
--        for Precision in reverse 6 .. 100 loop
--           Set_Default_Precision (Precision);
--           declare
--              A : MPFR_Float :=  To_MPFR_Float ("-1.8189894035458565e-12");
--              B : MPFR_Float :=  To_MPFR_Float ("2152019766.8571167");
--              C : MPFR_Float :=  To_MPFR_Float ("-914912179131808.0");
--              D : MPFR_Float :=  B * B - To_MPFR_Float (4) * A * C;
--              SD, NOM, DENOM, DIV, NDIV : MPFR_Float;
--           begin
--              Put_Line ("Precision:" & Precision'Img);
--              Put_Line ("D:" & Image (D));
--              SD := Sqrt (D);
--              Put_Line ("Sqrt(D):" & Image (SD));
--              NOM := B - SD;
--              Put_Line ("NOM:" & Image (NOM));
--              DENOM := (To_MPFR_Float (2) * A);
--              Put_Line ("DENOM:" & Image (DENOM));
--              DIV := NOM / DENOM;
--              Put_Line ("DIV:" & Image (DIV));
--              NDIV := To_MPFR_Float (0) - DIV;
--              Put_Line ("NDIV:" & Image (NDIV));
--
--              declare
--                 SD         : constant MPFR_Float := Sqrt (D);
--                 Solution_1 : constant MPFR_Float := -(B - SD) / (To_MPFR_Float (2) * A);
--                 Solution_2 : constant MPFR_Float := -(B + SD) / (To_MPFR_Float (2) * A);
--                 function Val (X : MPFR_Float) return MPFR_Float is
--                 begin
--                    return A * X * X + B * X + C;
--                 end Val;
--              begin
--                 Put_Line ("Solution_1:" & Image (Solution_1) & ", F(X)=" & Image (Val (Solution_1)));
--                 Put_Line ("Solution_2:" & Image (Solution_2) & ", F(X)=" & Image (Val (Solution_2)));
--              end;
--           end;
--        end loop;
--     end;

end MPFR_Example;
